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1. INTRODUCTION 

Development of new constitutive models for finite element applications 
represents a very important area of research in engineering disciplines. 

This is evidenced by research activities, for example, associated with high 
temperature composites [1,2], reinforced concrete [3], geotechnical materials 
[4,5], The efforts in constitutive research involve the development of 
mathematical relationships for predicting nonlinear response of materials, 
derivation of material stiffness matrix appropriate for finite element calcu- 
lations, computer implementation, and finally, coding verifications. Obvious- 
ly, the entire process requires a great deal of manual algebraic manipula- 
tions and computer programing. Hence, the response time for the related ef- 
forts is quite long, in the order of many months. As a result, it is rather 
difficult for the researcher to introduce any significant changes or modifi- 
cations into the constitutive theory, since the required effort is rather te- 
dious. Moreover, the outcome of research work 'may be affected by human errors 
which are often difficult to detect. In view of this discussion, it is ap- 
parent that symbolic manipulations can provide a significant incentive to the 
development of constitutive theories and their finite element applications. 

With the availability of MACSYMA or VAXIMA (i.e. VAX computer version of 
MACSYMA) , it becomes possible to derive the required material matrix of a 
constitutive model in an efficient way. The obvious advantages of using VAX- 
IMA are several: i) reduce manual tedium for deriving the material stiffness 
matrix, ii) improved reliability of analysis results, iii) quick response 
time for constitutive model development. However, direct application of VAX- 
IMA will not be trouble-free. One major obstacle is the exponential growth 
of algebraic expressions during intermediate derivations, which requires 
significant storage space and increased computer time. Moreover, it is also 
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possible to convert the derived expressions directly into Fortran coding. 

Then problems associated with modularity and interface with the main program 
must be addressed. 

Application of symbolic manipulations to finite element analysis is not 
new. Most of the previous work was concentrated, however, on the derivation 
of element stiffness and mass matrices [6-11]. No published work was found 
on the application of this procedure to the development of material matrices, 
although the general concept is somewhat similar. The objective of our re- 
search is to use symbolic manipulations for the derivation of a class of mat- 
erial matrices for finite element analysis; namely, elasto-plastic materials. 
The scope of our work includes derivation of material matrices and automatic 
Fortran code generation. In this paper, we will demonstrate a systematic ap- 
plication of the symbolic mathematical package, VAXIMA, the method of expres- 
sion simplifications, and code generation in the form of RATFOR. Three sam- 
ple constitutive models are included to illustrate the procedures developed. 
They are: von Mises metal plasticity, Drucker-Prager soil plasticity model, 

and a plasticity-based model for concrete. These models have been extensive- 
ly used for different finite element analyses in structural and geotechnical 
engineering fields. 
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2. THEORETICAL EQUATIONS 

For the sake of discussion, the stress-strain equations for elasto- 
plastic materials are briefly outlined in this section. More detailed des- 
criptions of these relations can be found in several texts [e.g., 12-14]. It 
is noted that our primary purpose here is to derive the general form of the 
material stiffness matrix as commonly used in the displacement-based finite 
element analysis. 

The first basic assumptions in the incremental (flow) theory of plasti- 
city is the additive decomposition of the total incremental strain vector, 

e p 

de, into elastic and plastic components, de and de , respectively. In ad- 
dition, the incremental elastic strain components are often assumed to be 
linearly related to the increment of stress vector (generalized form of 
Hooke's Law), 

da = C (de - de^) (1) 

#v /v /v 

where is an elastic material stiffness matrix. 

Thus, the main task in the formulation of the elasto-plastic model is 
concerned with establishing the manner in which the plastic strain increments 
are related to the stress increment vector and the history of deformation. 

To this end, three fundamental assumptions of plasticity theory are employed. 
These are: i) the yield (loading) function defining the limit of elasticity 

of the material during the course of plastic deformations ii) an appropriate 
hardening rule specifying the manner for the evolution of the yield surfaces 
during plastic straining (e.g. isotropic, kinematic, or mixed hardening, etc) 
and iii) a flow rule that provides the general form of the incremental plas- 
tic stress-strain relationships (e.g. associated flow or normality rule, or 
the non-associated flow rule). 
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Under isothermal conditions, the yield function is expressed as 

f - f(a , k) (2) 

where a is the stress vector, and k represents a strain-hardening parameter 
which may vary as a function of plastic deformation history or other state 
variables. Note that, in general, one or more strain-hardening parameters 
may be used, and these may actually be scalars {e.g., accumulated plastic 
work) or tensorial quantities (e.g., plastic strain components). However, 
for convenience, we only use one scalar hardening parameter here, since all 
of the specific plasticity models considered are of this type. Also, asso- 
ciated flow rule is employed in the three models discussed. 

Adopting the normality rule, the plastic flow (or increment of plastic 
strain vector) is given by 

de p * dX || (3) 


where dx is a positive scalar quantity often referred to as the loading para- 
meter or plastic multiplier, which generally depends on the current state of 
stress cr, incremental stresses do, and loading history. 

Based on the above relationships, and employing the so-called consis- 
tency condition [13], one can easily derive the general form of the incre- 
mental stress-strain equations for a material undergoing plastic deformations 
[12,14]; that is 


where C designates the familiar elasto-plastic material stiffness matrix 
which has the form 



P 

where C is a plastic matrix given by 

'V 


( 5 ) 
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In addition to the relationship in Eq. (4), the matrix c" is also used for 

T E P 

the evaluation of element stiffness matrix jc = /// B C B dv. 

It is seen from Eq. (6) that in order to obtain specific expression for 

the elasto-plastic matrix C ^ one has to manipulate the derivatives of 
yield function with respect to a and * and then carry out appropriate matrix 
multiplications. For a complex mathematical expression of f, the associated 
manipulation can be quite tedius if it is done manually. In the next sec- 
tion, we outline the procedure through which this can be done conveniently by 
symbolic manipulations using VAXIMA. 


3. SYMBOLIC MANIPULATIONS 

In order to find the explicit expression of elasto-plastic matrix, i.e., 
Eq.(5), for a given material model, four types of mathematical manipulations 
need to be employed; i) differentiations of the yield function with respect 
to the stress or other variables, ii) matrix multiplications, iii) grouping 
of like-terms, and iv) expression simplifications. It was pointed out ear- 
lier that in most cases the results obtained from direct application of VAX- 
IMA would not be useful due to the problem associated with expression growth. 

For this reason, a strategy must be developed to obtain an optimal (or simp- 

EP 

lified) form for the matrix C . The essence of our strategy consists of: 

1. A structured derivation procedure to avoid redundant manipulations 
and to minimize expression growth. 

2. Factorization and expression simplification through user interven- 
tion with interactive coding. 

3. Introduction of intermediate variables. 

4. Taking advantage of permutation and symmetry relationships of the 
terms and matrix involved during intermediate derivations. 

With the above guidelines in mind, the derivation of elasto-plastic 
material matrix, i.e. Eq. (5), involves the following: 

1. Finding the derivatives 

JB 

and 

q 


/ 8f 3_f 3f 

3o ll* 9o 22 ’ ^°33 


3f 3f 3f » 
’ 3 o 12* 3o 23’ 3o 31 
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p 

2. Performing matrix multiplications for the numerator of C : 


E 

C • P 


and for the denominator: 

T E T 

P C P and p • q 

IM (V 

3. Conducting expression simplifications during the course of derivations. 
Two simplification condi tons are often used for substitutions: 


S 11 + s 22 + S33 = 0 


and 


2 2 2 2 2 2 , 

S + S + S + 2(S + S + S ) = 2J 

11 22 33 12 23 31 2 


(7) 

(8) 


In the sequel, three specific examples of plasticity material models are 
employed to demonstrate our procedure outlined ( in the above. 


3.1 von Mises Metal Plasticity Model 

We consider first the von Mises model with isotropic strain hardening 
for metal plasticity as an example to demonstrate our procedure. In this 
case, the yield function f is given by 

f = j § T S - k 2 (9) 

where S represents the vector of stress deviators and 

S = a - a • 6 (10) 

~ m ~ 

6 is a vector of Kroneckle delta 

6 = ( 1 , 1 , 1 , 0 , 0,0 ) ( 11 ) 
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( 12 ) 


and a is the mean stress, a scalar quantity given by 
m 

°m * T (a ll + a 22 + a 33 * 


The parameter k is a function of plastic work 


k = k(W p ) 


Also note that k is related to the effective stress a by [12] 

e 

_ 1 _ 


(13) 


(14) 


Using VAXIMA, at first we evaluate 


and 


p = (S ,S ,S ,2S ,2S ,2S ) 

~ 11 22 33 12 23 31 


d = H(o ,a ,o ,a ,a ,a ) 
11 22 33 12 23 31 


(15) 


(16) 


In the above, the common factor H is the so-called plastic modulus and it 
is found to be 


2 3f 

7 a e 3^ 


(17) 


Next, evaluating 


Y = C 


* ~ n+v)(i-2v) 


(l-v)Sn+vS22 + vS33 
vSn+(l-v)S22+vS33 
vSi i+vS9?+(l-v)S33 
(l-2v)S 12 
(1-2v)S23 
(1-2v)S 3 i 


(18) 


The numerator of C is equal to 


T 2 

Y • Y = 4G • S 

~ ~ -2 


(19) 
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where G is the shear modulus, and S 2 is a 6 x 6 matrix in which the entries 
are the products of stress deviators, i.e.. 


S 2 = S • S T (20) 

rv ft# fV 

P 

Then, we evaluate the denominator of C . For this purpose, we utilize the 

ft# 

two simplification conditions given in Eqs. (7) and (8). Thus, the follow- 
ing simplified expressions can be found 


and 


T E 
p C p 

ry# 




2 

1 



( 21 ) 

( 22 ) 


By combining Eqs (19), (21) and (22) with Eq. (5), we finally obtain 


.EP 


C E - 


3G 


Og (H + 2G) 


(23) 


The above expression corresponds to that given in [5]. 


3.2 Orucker-Praqer Soil Plasticity Model 

We now consider a more complex material model; namely, the Drucker- 
Prager, perfect-plasticity model used extensively for geotechnical materials 
In addition to the above procedures, intermediate variables have to be intro 
duced in this case in order to avoid the problem of expression growth. The 
yield function of Drucker-Prager model assumes the form [16] 

1/2 

f-j +ol-k (24) 

2 1 

where 1^ is the first stress invarient; J 2 is the second invarient of stress 

deviators; a and k are material constants. 
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If we follow the same procedure for the von Mises model without intro- 
ducing any intermediate variables, growth of algebraic expressions becomes 

P 

apparent. For example, the first three entries in the first row of C are 

listed in Fig. 1. There are twenty-one similar entries when the symmetry 
P 

condition of C is taken into consideration. However, after we have 
introduced the intermediate variables a and ff , 

a = (ai, a 2 .... ag) (25) 

ff * (ffi, ff 2 ,...ffe) (26) 

where a^, ff^ are defined in Fig. 2, then 

EP T 

C = ff • (ff) (27) 

fV IV»V AAV 

Of course, in computer coding we only need to perform matrix multiplication 

EP 

for either the upper or lower half of C owing to its symmetry properties. 

3.3 Concrete Plasticity Model 

The final illustrative example to be considered here is a concrete plas- 
ticity model proposed recently by Chen and Chen [13,17]. The derivation of 
elasto-plastic matrix for this hardening plasticity model is quite tedius due 
to its complex mathematical expressions. In fact, through the use of symbol- 
ic manipulations an error was found in the coefficient of the plastic matrix 
published in the literature in [12,13,17]. The error is associated with 
the shear stress terms, which usually do not appear for simple test cases 
like simple compression or biaxial compression tests that have been used in 
various model verifications. Moreover, the error terms do not appear when the 
model is reduced to the von Mises theory. 

The yield function of the concrete plasticity model is given by 

10 . 


f 


(28) 


J 2 ♦ 


t h + " 


1 " 7 l l 


x 2 = 0 


where a and B are material constants; t is an effective stress; n, a parame- 
ter, whose definition varies with the stress state: n = 0 for compression - 
compression stress states; n * -1/6 for tension - compression or tension - 
tension stress states. 

For this model, the problem of expression growth becomes prohibitive if 
intermediate variables are not introduced. By successive manipulations with 
VAXIMA the following variables are identified: 


a = ^ ( p • « + & 

( 29 ) 

S = (S U , S22, S33, 2Si2, 2S23, 2S31) 

( 30 ) 

5 = (1, 1, 1, 0, 0, 0) 

( 31 ) 

w ' p l« * *1 

( 32 ) 

C 


u = £ C 2 J 2 * 30 2 + 2<sj * S 2 2 * S 3 2 )] 1/2 

( 33 ) 


v = -^ [2J 2 + 3p 2 + v( 4J 2 + 3p 2 )] (34) 

m 

E c = (1W)(1-Tv7 {35) 
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(36) 


P * nil + y (B + ax 2 ) 

m = 1 - y I| (37) 

where E = Young's modulus, v * Poisson's ratio, and H = plastic work harden- 
ing paratmeter. It is noted in the above that the underlined terms for u 
were missing in the published expression [12]. The addition of these terms 
was verified by both symbolic manipulations using VAXIMA and independent man- 
ual derivations. 

Moreover, we introduce a vector ff(i), i=l, 2, ...» 6, of which the 
first three components are given by 

ff ( i ) * (I + v ft) a(i), i=l, 2, and 3 

and the last three components are 

ff (i ) = i 1 Z — ) ■ a(i ) , i=4, 5 and 6. 

IV~ C /V 

where I is a 3 x 3 identity matrix and 



The above procedure has been written in LISP program language with di- 
rect access to the internal data structure of VAXIMA. Hence, the package can 
be used for the symbolic manipulations of any elasto-plastic material matrix 
with two special features: i) no expression growth problem, and ii) high ef- 
ficiency in terms of CPU time. 
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4. AUTOMATIC CODE GENERATION 

There is a definite advantage to convert the generated symbolic expres- 
sions into FORTRAN statements for finite element computations. By doing 
so, not only can the manual effort be avoided, but also it provides in- 
creased reliability on the constitutive relations. Instead of generating 
FORTRAN directly, we have utilized a generator called GENTRAN (symbolic to 
numerical code GENerator/TRANslator)[18] which has the ability to produce a 
RATFOR or C program in the form of a subroutine or part of a subroutine. 
Subsequently, the FORTRAN statements are generated from the RATFOR through a 
preprocessor. 

Several systems are available in converting symbolic expressions to FOR- 
TRAN statements, such as MACTRAN [19], VAXTRAN [20] and REDUCE [21]. The 
MACTRAN Package converts MACSYMA equations and other expressions into FORTRAN 
code, and provides a text processor which allows the derived FORTRAN code 
segments to be interspersed with fixed code from program skeletons. Similar 
features are given in the REDUCE and VAXTRAN systems, except that VAXTRAN 
was written specifically for VAXIMA. All these packages represent a first 
step towards providing an interface between symbolic manipulations and nu- 
merical computations. However, they do not provide a convenient way to gen- 
erate statements such as declarations, control-flow structures, I/O state- 
ments, functions, and subroutines. These statements, in general, are neces- 
sary for generating a complete and efficient FORTRAN program. For this rea- 
son, we have chosen to use a package called GENTRAN which was written in 
FRANZ LISP under the VAXIMA environment. 

The immediate concern in generating a subprogram to interact with a 
finite element code is the interface problem. To minimize such problems, we 
have designed a template file shown in Fig. 3 which is somewhat universal for 
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various finite element codes. For a specific plasticity material model (such 
as the von Mises or Drucker-Prager model), the material matrix subroutine is 
completed by including the generated statements in the template file as indi- 
cated in Fig. 3. 

To demonstrate how the RATFOR code is generated by GENTRAN, we consi- 
der again the von Mises model. First, we define the S 2 matrix according to 
Eq.(18). Secondly, the elasto-plastic matrix is evaluated from Eq. (23). 

Let 

S(i) = components of stress deviator, i =1, 2, ...» 6. 
CE(i,j) = elastic material matrix 
CEP ( i , j ) = elasto-plastic matrix 
CP(i,j) = plastic material matrix 

FACTOR = -M 

c/(H+2G) 

Then the RATFOR code for the von Mises model is given as follows: 


1. 

for ( i =1 ; 

i<=6; 

i=i+l) 

2. 

for ( j=i ; 

j<=6; 

j=j + l) 

3. 

CP(i,j) = 

FACTOR 

* S(i) 

4. 

for (i=l; 

i<=6; 

i=i+l) 

5. 

O 

"5 

Cj. 

II 

V • 

j<=6; 

j=j+l) 

6. 

CEP(i.j) 

» CE(i , j 

) - CP( 

7. 

CEP(j.i) 

= CEP(i, 

j) 


The translated FORTRAN code can be found in the Appendix. 

The RATFOR code of the Drucker-Prager model is slightly different from 
that of the von Mises model due to the use of intermediate variables. In 
this case, let 
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FF(i) - intermediate variables as defined in Eq.(26), i*l,2,...,6 
FACTOR = 1/ (G + «B) 

Bs l“ 6 (rar) 

V = Poisson's ratio 
Then the RATFOR code is 

1. for (1*1; i<=6; i=i+l) 

2 - f°r ( j=i ; j<=6; j=j+l) 

3 - CP(i,j) = FACTOR * FF(i ) * FF(j) 

4. for (1=1; i<=6; 1=1+1 ) 

5 - for (j=i; j<=6; j=j+l) 

6 - CEP(i.j) = CE(i.j) - CP(i,j) 

7. CEP(j,i) = CEP(i,j) 

The translated FORTRAN code for the Orucker-Prager model is given in the 
Appendix. 

Finally, with the introduction of the intermediate variables in Eqs. 

(29) (37), which were obtained through the factorization of VAXIMA, the 
RATFOR coding of the concrete model becomes identical to that of the Drucker- 
Prager model. The corresponding FORTRAN code is listed in the Appendix. 


15 



5. CONCLUSION 


A systematic procedure to perform symbolic manipulations using VAXIMA 
and FORTRAN code generation of elasto-plastic material matrices for finite 
element applications has been developed. The unique features of the proposed 
procedure are: i) the problem of expression growth was alleviated by intro- 
ducing intermediate variables and step-wise expression simplifications, ii) 
the material matrix is automatically converted into FORTRAN coding, and iii) 
the use of a template file to ease the interface problem. This procedure can 
be applied not only to plasticity models with associated flow rules, but also 
to models with non-associated flow rules. 

The potential benefits of the proposed procedure are two-fold: i) it 

can avoid manual tedium for constitutive model development, and ii) it pro- 
vides increased reliability on the model for finite element applications. 

The same concept can be extended to other types of constitutive model devel- 
opment. For example, in the finite element analysis of viscoplastic consti- 
tutive models, the formation of Jacobian matrix for numerical integration 
requires lengthy algebraic manipulations of the rate stress-strain equations. 
Such manipulations can be easily performed by a well-designed VAXIMA proce- 
dure. Once the mathematical relations are derived, automatic code generation 
should become apparent. 
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C 11 = ^ 4a2j 2 ‘ ^ J 2 S 22 + S 22 ~ 4ot J 2 S 33 


2 2 2 
+25 S + S + 8a v J + 4av J S nn - 4v S 
22 33 33 2 2 22 22 


2 2 2 

+ 4av J S - 8v S„ - 4v S„ + 4a v J 
2 33 22 33 33 2 


2 2 2 2 
+ 8av J S + 4v S + 8 av ■ J , S 

2 22 22 2 33 


2 2 2 
+ 8 v S S +4vS ] 
22 33 33 


C 12 = 4J 2 ^ 4a J 2 ’ S 22 ‘ 2a J 2 S 33 “ S 22 S 33 


2 2 

+ 8a v J + 4vS + 2a v J S + 4vS S 
2 22 2 33 22 33 


2 2 2 2 2 2 
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2 2 2 2 2 2 
+ 4avJ -4v S + 4a v J S - 4v S S 

2 33 2 22 22 33 


etc. 

Figure 1 Typical Entries of Plastic Matrix of Drucker-Praeger Model 
Without the Use of Intermediate Variables. 
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a l = “ + “2^" S 11 

a 2 = ° + S 22 

a 3 = ° + S 33 

- 1 c 

a 4 b 12 

, _ 1 c 

a 5 " b Z3 

, _ 1 C 

a 6 " J 2 b 31 

ff = (l-2v) a + v a , i = 1, 2, or 3 
i i o 

f f . = . aj , j = 4, 5, or 6 


a 

o 


a 

1 


+ a 

2 


+ a 
. 3 


Figure 2 Intermediate Variables for Expression Simplification of 

Drucker-Praeger Model 
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SUBROUTINE EPMTRIX (STRESS, STRAIN, CEP, H) 

****************************** 

* * 

* A TEMPLATE FOR ELASTO-PLASTIC MATERIAL MATRIX * 

* * 
****************************** 

CEP - ELASTO-PLASTIC MATRIX; H - STRAIN HARDENING PARAMETER 

E - YOUNG'S MODULUS ; PV - POISSON'S RATIO 

CONI , C0N2, ADDITIONAL MATERIAL CONSTANTS 

COMMON /MDATA/ E, PV, CONI , C0N2, C0N3, CON4, CONS 

DIMENSION STRESS (i> , STRAIN < i > , CEP (6, 6) , CE <6, 6) ,CP(6,6> 

DIMENSION A (6) , FF (6) 

DEFINE ELASTIC MATERIAL MATRIX 


DO 10 I 
DO 10 J 
10 CE ( I , J ) 
COEF 
CE ( 1 , 1 > 
CE (1,2) 
CE ( 1 , 3 > 
CE (2,2) 
CE (2,3) 
CE (3,3) 
CE (4,4) 
CE (5,5) 
CE (6,6) 

DO 20 I 
DO 20 J 
20 CE(J, I) 


= 1,6 
= 1,6 

= O. 

= E/( (1. + 
= COEF* ( 1 . 
= COEF*PV 
= COEF*PV 
= CE ( 1 , 1 ) 

= CE (1,3) 

= CE ( 1 , 1 ) 

= COEF* ( 1 . 
= CE (4,4) 

= CE (4,4) 

= 1,6 
= 1,6 
= CE ( I , J ) 


PV>*(1. - 2. *PV) ) 
- PV) 


- 2. *PV)/2.0 


DEFINE STRESS DEVIATORS 


SIGM = (STRESS (1 ) +STRESS ( 2 ) +STRESS ( 3 ) >/3.0 

511 = STRESS ( 1 ) - SIGM 

522 = STRESS (2) - SIGM 
S33 = STRESS (3) - SIGM 

512 = STRESS (4) 

523 = STRESS (5) 

S31 = STRESS (6) 


( ) 

( GENERATED FORTRAN CODE ) 

( ) 

RETURN 

END 

Figure 3 A Template File for Elasto-Plastic Material Matrix 





APPENDIX 


C 

C VON MISES METAL PLASTICITY 

C 

C DEFINE MATERIAL PARAMETERS 

C 

G =E/(2.0*(1+PV) ) 

XJ2 = (SU*SU+S22*S22+S33*S33+2*(S12*S12+S23*S23+S31*S31> >/2.0 
FAC =G/ (XJ2* (H+2*G) > 

A ( 1 ) =S 1 1 
A <2) =S22 
A (3) =S33 
A <4) =S12 
A (5) =S23 
A (6) =S31 
1 = 1 

2300 IF ( . NOT. ( I . LE. 6) ) GOTO 2302 
J=I 

2303 IF (.NOT. (J.LE.6) )GOTO 2305 
CP(I, J)=FAC*A(I) *A(J) 

J=J+1 
GOTO 2303 

2305 CONTINUE 
1=1 + 1 
GOTO 2300 

2302 CONTINUE 
1 = 1 

2306 IF (.NOT. (I.LE.5) >GOTO 2308 
J=I + 1 

2309 IF (.NOT. (J.LE.6) ) GOTO 2311 
CP(J, I)=CP(I, J) 

J=J+1 
GOTO 2309 

2311 CONTINUE 
1 = 1 + 1 
GOTO 2306 

2308 CONTINUE 
1 = 1 

2312 IF (.NOT. (I. LE. 6) ) GOTO 2314 
J=1 

2315 IF(. NOT. (J.LE.6) ) GOTO 2317 
CEP ( I , J ) =CE ( I , J ) -CP ( I , J ) 

J=J+1 
GOTO 2315 

2317 CONTINUE 
1 = 1 + 1 
GOTO 2312 

2314 CONTINUE 



nnnnn 


DRUCKER-PRAGER MODEL 
DEFINE MATERIAL PARAMETERS 
AFA =CONl 

COEF=E/ ( (l+PV) * (1-2*PV> ) 

XJ2 = (S11*S11+S22*S22+S33*S33+2*(S12*S12+S23*S23+S31*S31> >/2.0 
PQ =0 

QTCQ=C0EF/2* ( ( 1-2*PV> +6*AFA**2* ( l+PV) ) 

WINV=QTCQ+PQ 
W =C0EF**2/WINV 
A ( 1 ) =AFA+1 /SQRT (X J2> *S1 1 /2. O 
A < 2) =AFA+ 1 /SORT < X J2 > *S22/2 . 0 
A (3> =AFA+1 /SORT (X J2) *S33/2. O 
A (4) = 1 /SQRT < XJ2) *S12 
A(5)=1/SQRT(XJ2) *S23 
A(6>=1/SQRT(XJ2) *S31 
SUM = A ( 1 ) + A < 2 > + A ( 3 ) 

FF ( 1 ) = < 1 -2*PV> *A (1 ) +SUM*PV 
FF (2) = ( 1— 2*PV) *A(2)+SUM*PV 
FF (3> = < i— 2*PV) *A(3)+SUM*PV 
FF <4) = ( 1— 2*PV) *A<4> /2.0 
FF(5)=(1-2*PV) *A(5)/2.0 
FF (6) = ( 1-2#PV) *A(6)/2.0 
1 = 1 

2300 IF (.NOT. (I.LE.6) ) GOTO 2302 
J=I 

2303 IF (.NOT. (J.LE.6) > GOTO 2305 
CP ( I , J ) =W*FF ( I ) *FF ( J > 

J=J+1 
GOTO 2303 

2305 CONTINUE 
1 = 1 + 1 
GOTO 2300 

2302 CONTINUE 
1 = 1 

2306 IF (.NOT. ( I . LE. 5) > GOTO 2308 
J=I + 1 

2309 IF (.NOT. (J.LE.6) ) GOTO 2311 
CP ( J , I ) =CP ( I , J ) 

J=J+1 
GOTO 2309 

2311 CONTINUE 
1 = 1 + 1 
GOTO 2306 

2308 CONTINUE 
1 = 1 

2312 IF (.NOT. (I.LE.6) ) GOTO 2314 


A-2 



J=1 

2315 IF (.NOT. (J.LE.6) )GOTO 2317 
CEP ( I , J ) =CE < I , J ) -CP ( I , J ) 
J=J+1 
GOTO 2315 
2317 CONTINUE 
1 = 1 + 1 
GOTO 2312 
2314 CONTINUE 


A-3 
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CONCRETE PLASTICITY MODEL 

DEFINE MATERIAL PARAMETERS 

XA =CONl 
XB =CON2 
XN =C0N3 
XM =1-XA*SIGM 
COEF=E/ ( < 1+PV) * < 1-2*PV> ) 

SIG3=3*SIGM 

XJ2 * (SI 1*S1 l+S22*S22+S33*S33+2* (S12*S12+S23*S23+S31 *S31 ) ) /2. 0 
TU =SQRT( (XJ2+XN*SIGM*SIGM/2.0+XA*SIGM/3.0)/XM> 

HP =2*TU*H 

RO =XM*SIGM+(XB+XA*TU*TU>/3.0 

U =(HP/XM)*SQRT(2*XJ2+3*R0**2-2*(S12**2+S23**2+S31t*2> ) 

V =COEF/ (XM*XM> * (2*XJ2+3*P**2> * ( 1-2*PV> 

WINV=U+V 

W =C0EF**2/WINV 
A ( 1 ) = (RO+S1 1 > /XM 
A ( 2 > = ( R0+S22 ) / XM 
A (3) = (R0+S33) /XM 
A<4)=2*S12/XM 
A (5) =2*S23/XM 
A (6) =2*S31/XM 
SUM = A ( 1 ) +A ( 2 ) +A < 3 ) 

FF ( 1 ) = ( 1-2*PV> *A ( 1 > +SUM*PV 
FF (2) = ( 1-2*PV> *A < 1 ) +SUM*PV 
FF <3) = ( 1-2*PV> *A (3) +SUM*PV 
FF(4)=(l-2*PV)*A(4>/2.0 
FF (5) = ( 1-2*PV> *A (5) /2. O 
FF<6)=(1-2*PV> *A(6>/2.0 
1 = 1 

2300 IF < . NOT. ( I . LE. 6) > GOTO 2302 
J=I 

2303 IF ( . NOT. ( J . LE. 6) ) GOTO 2305 
CP ( I , J ) =W*FF ( I ) *FF ( J ) 

J=J+1 
GOTO 2303 

2305 CONTINUE 
1=1 + 1 
GOTO 2300 

2302 CONTINUE 
1 = 1 

2306 IF (.NOT. (I. LE. 5) ) GOTO 2308 
J=I + 1 

2309 IF ( . NOT. ( J. LE. 6) ) GOTO 2311 
CP(J, I)=CP(I, J) 

J=J+1 


A-4 



GOTO 2309 

2311 CONTINUE 
1 = 1+1 
GOTO 2306 

2308 CONTINUE 
1=1 

2312 IF ( . NOT . < I . LE. 6) ) GOTO 2314 
J=1 

2315 IF < . NOT. ( J . LE. 6) ) GOTO 2317 
CEP ( I , J > =CE ( I , J ) -CP < I , J > 
J=J+1 
GOTO 2315 

2317 CONTINUE 
1 = 1 + 1 
GOTO 2312 

2314 CONTINUE 
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